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ABSTRACT 

Mean-motion resonances (MMRs) are likely to play an important role both during and after the 
lifetime of a protostellar gas disk. We study the dynamical evolution and stability of planetary systems 
containing two giant planets on circular orbits near a 2:1 resonance and closer. We find that by having 
the outer planet migrate inward, the two planets can capture into either the 2:1, 5:3, or 3:2 MMR. 
We use direct numerical integrations of 1000 systems in which the planets are initially locked into 
one of these resonances and allowed to evolve for up to ~ 10^ yr. We find that the final eccentricity 
distribution in systems which ultimately become unstable gives a good fit to observed exoplanets. 
Next, we integrate ~ 500 two-planet systems in which the outer planet is driven to continuously 
migrate inward, resonantly capturing the inner; the systems are evolved until either instability sets in 
or the planets reach the star. We find that although the 5:3 resonance rapidly becomes unstable under 
migration, the 2:1 and 3:2 are very stable. Thus the lack of observed exoplanets in resonances closer 
than 2:1, if it continues to hold up, may be a primordial signature of the planet formation process. 
Subject headings: celestial mechanics — planetary systems: formation — planetary systems: proto- 
planetary disks — planets and satellites: general 



1. INTRODUCTION 

The discovery of extr asolar planets around su n-like 
stars (|Mavor fc Queloz I 119951 : iMarcv fc Butler I Il995l 
Il998f) has revealed that early large-scale migration likely 
plays an important role in the formation of planetary 
systems. For example, observations have revealed a 
nontrivial number (roughly 6% overall) of close orbit- 
ing gas giants, called 'hot Jupiters,' planets that would 
have great difficulty formin g in their present locations 
(jBodenheimer et al. I |200C1[) . In fact, disk-planet in- 
teraction theory predicts disturbingly short migration 
time scales that seem to threaten the survival of ev- 
erything ran ging from planetary embryos to gas giants 
([Ward II1997I ). Additionally, there are currently at least 
eight planetary sys tems that have two planets in MMR 
(lUdrv et al. I (20071). Thes e include syste ms such as 



GJ 876 (IMarcv et al. l[200lb and HD 82943 (iLee et al. I 
I2006D with two planet s in a 2:1 resonance, and HD12661 
([Fischer et al. 1 120031 ). which contains two planets in 
what may be a 6:1 resonance. The origins of such reso- 
nant systems are thought to be due to disk-planet inter- 
actions that induce angular momentum transfer and dif- 
feren t ial migration (ISnellgrove et al. "20011: ILee fc Peale I 
I2002t iPapaloizou I l2003: Kley et al. 2004). 

There are currently two competing theories to explain 
planet formation. The "core accretion" model starts 
with the sedimentation and coUisional growth of dust 
grains a nd smaller plan etesimals in the protoplanetary 
disk (see iLissauer |[l993l and references therein). These 
rocky cores (protoplanets) continue to build up until 
their gravitational influence allows for the accretion of 
the surrounding gas, forming gas giants. The other the- 
ory i nvokes direct gravitational i nstabilities in the disk 
(^e.g.. [Cm^ero^[l97llBo^ 



In these theories the dynamical relaxation time of 
a planetary system can, in principle, be longer than 
the time scale for planet formation. Therefore, long 
term stability is not guaranteed for a planetary sys- 
tem's initial configuration. Gravitational interactions 
between massive protoplanets can lead to orbit cross- 
ing and instability, resulting in m assive planets being 
thrown closer to their parent star s (|Rasio fc Ford ifigga 
IWeidenschilling fc MarzaiT1ll996l ). In the case of one of 
the planets being ejected, the surviving planet is left with 
high eccentricity (some with e > 0.90). Additionally, 
migration occurs within the disk due to the exchange 
of angular momentum between the planet and the sur- 
rounding disk material through planet-disk interactions 



(iGoldreich fc TremaTnel I197L. 
119791 119931 : IPapaloizou fcLhTl 
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'Lin fc Papa loizou I 
Ward 1986). In a 
laminar disk, there are two basic modes of migration: For 
an embedded body, imbalance between the torques from 
the inner and outer parts of the disk is thought to lead to 
orbital decay ' this is comm only referred to as type I mi- 
gration fe.g.. lWard 1119971 ). If a body IS massive enougn, 
it locks itself to the disk by opening a deep annular gap, 
and is thus carried along as the disk ac cretes onto the star 
(|Lin fc Papaloizou lll986l : IWard Ill997f ): this is called type 
II migration.^ When orbital migration occurs for differ- 
ent planets at different rates, convergent migration and 
locking into mean-motion commensurabilities can occur. 

For the two planet case, the regions of dynamical 
stability (i.e., the region where close encounters are 
impossible) can be determined analytically by study- 
ing how the value of the Jacobi integral relates to 
the t opological (Hill) stability of the system (jGladman I 
Il993f ). However, this criterion does not take into ac- 
count the fact that planets may be in resonance. It 
has been shown that resonant systems have additional 

There is also a proposed type III migration, which involves 
disk material flowing through the co- orbital region of the pla net, 
which will not be considered here (see IPapaloizou et al. ir2007D . 
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regions of stability (|Barnes fc Greenberg Il2007[) . which 
current analytic criteria pre dict to be unstable. In fact, 
iBarnes fc Greenberg I ()2007f ) have shown that nearly all 
observed resonant systems lie in these extended regions. 

Overall, planet-planet-disk interactions are likely key 
in determining the final configuration of a planetary 
system. To date, each type of dynamical process 
(planet-disk interactions, dynamical instabilities, and 
resonances) has usually been considered individually. 
The motivation of this study is to explore physical situ- 
ations where all of these processes occur simultaneously. 
In particular, we consider gap-opening planets located 
close to particular resonances within a protoplanetary 
disk. This results in convergent migration and reso- 
nance capture, but possibly also eventual instability due 
to close encounters. 

There are two objectives to this study. The first is to 
better quantify the regions of stability for particular res- 
onances. More explicitly, we will consider the 3:2, 5:3, 
and 2:1 resonance for two gap-opening planets of varying 
mass. We place an upper bound on the possible masses 
for planets in these three resonances. This provides an 
interesting problem in dynamics by allowing us to more 
generally map the parameter space for these particular 
resonanc es in and out of the p rotoplanetary disk, com- 
pared to iBarnes fc Greenberg I ([2007) who studied par- 
ticular observed extrasolar systems. Additionally, in the 
case where there are instabilities, this study allows us 
to see how the final configurations (e.g., distribution of 
eccentricity) differ, if at all, fro m non-rcsonant planet 
scattering (see, e.g.. iFord et al. I 2001: Chattcricc ct al. 
[200l . 

Our second objective is to study the dynamics of gas 
giants in the region of the protoplanetary disk where we 
believe these planets form, between 5 and 10 AU. The 
three above resonances are typical resonances that plan- 
ets can resonantly capture into, since they require little 
to no eccentricity to give a high probability of capture 
when planets migrate at rates consistent with Type II 
migration. Here we are assuming that planets are born 
on nearly circular orbits. We then study the evolution 
as the coupled system of two planets migrates with the 
disk. To date, the closest resonance to be observed — that 
is, the resonance with the smallest ratio of outer to inner 
period — is the 2:1. We aim to better quantify why we 
have not found any systems in the 3:2 or 5:3 resonances 
and measure how common we should expect each of these 
resonances to be in extrasolar planetary systems. 

This paper is organized as follows. In §2, we explain 
the numerical treatment used in both the resonance sta- 
bility and disk simulations. In §3, we give further details 
of the theory used in our models for studying the stabil- 
ity of the aforementioned resonances. We then present 
the results of those simulations. In §4, we develop the 
theory used to model the planet-disk interactions in the 
protoplanetary disk as well as the results of those runs. 
We conclude with a summary and discussion in §5. 

2. NUMERICAL METHODS 

Newton's equations of motion are integrated directly 
in the barycen tric frame using a f ifth-order Runge-Kutta 
(R-K) scheme ()Press et al. 1120 071. Although these direct 
integrations are computationally more taxing compared 
to symplectic integrations, they are able to handle ar- 



bitarily short dynamical times (e.g., when the planets 
are orbiting within a few solar radii of the star). Since 
our integrations run at most a length of 10^ yr, we al- 
low the R-K scheme to incur a maximum error of order 
1 part in 10^ per time step, where the time step in the 
three-body problem is variable. This ensures that the 
total accumulated error in energy, characterized through 
the ratio /S.E/Ei, where Ei is the initial system energy 
and Ai? is the difference between this value and the fi- 
nal energy of the system, remains less than 10~^ in the 
conservative case. Similarly, llA/Ai remains less than 
5 • 10~^, where Ai is the system's initial angular momen- 
tum, ensuring that the accumulated errors do not affect 
the final result. 

In all our simulations, we consider a two-dimensional 
three-body system, consisting of a solar-mass star and 
two planets with masses ranging between 0.25 Jupiter- 
masses (Mj) to 12 Mj. The system evolves until either 
a planet is ejected, a collision between any two bodies 
occurs, or the system evolves for 10^ yr. An ejection is 
defined as a planet being at least 250 AU from the star 
with positive energy. A collision between the star and 
a planet occurs when the relative distance between the 
two is less than 1.1 Rq, where Rq is the radius of a 
solar-mass star.'* Similarly, a collision between two plan- 
ets occurs when the relative distance is less than 2 Rj. 
The simulation continues for several years after the col- 
lision or ejection before ending the run. In the event of 
a collision, the pair is replaced with a single point mass 
with the same linear momentum as the original pair and 
a mass as simply the sum of the original masses. 

In these simulations, the units are selected so that 1 AU 
= 1, Mq = 1, and G = 1, where G is the gravitational 
constant. The planet's initial orbits are always set to be 
circular. 

3. RESONANCE STABILITY 

As mentioned in Section 1, the region of dynamical 
stability is sharply defined for t he two planet ca se, and 
its location is analytically known (jGladman 1119931 ). That 
is, given the masses of the two planets, the minimum 
orbital separation in the semi-major axes to prevent close 
encounters can be determined. Since we wish to consider 
two planets in a particular resonance, we can reverse the 
criterion to give analytically determined limits on the 
masses of the planets. We compute the separation in 
semi-major axes via Kepler's third law. Since the ratio 
of the orbital periods To/Ti > 1, where the subscripts 
denote the outer and inner planet, respectively, we get 
the following relation: = {To/Ti)^/'^ai. Therefore, the 
orbital separation normalized by the inner planet's semi- 
major axis is simply 

where we have defined p = To/Ti . So A is a constant for a 
given resonance. Now the analytic criterion for stability 
for low-eccentricity planets is given by 

A„„ = 2 • 31/6(^1 + /Z2)'/' + 2 • 31/3(^1 -I- ^2)2/3 

11^1 -I- III2 

3"/6(;,l+;,2)l/3+---' 

1 Rq + 1 Rj ~ 1.1 Rq, where Rj is the radius of Jupiter. 
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where fii ( = Mi/M(t,) is ch osen to be the larger of the 
two bodies (|Gladman II1993D . By setting A = Ami„, we 
can derive a maximum hmit on the masses of the planets 
to ensure stability, as suggested by Equation [21 Table 1 
gives some results for equal mass planets. Additionally, 
it compares the results whether one uses only one, two or 
three terms in the expansion (Equation [2]). One can see 
for resonances with larger orbital separations, truncating 
the series with the first term is inadequate to accurately 
describe the analytical limit. In fact, the addition of 
the second term results with over 100% correction to the 
mass limit for resonances above the 2:1 MMR. 

3.1. Initial Conditions 

We randomize the initial semi-major axis of the inner 
planet between 5 and 6 AU. For the 3:2 resonance, the 
masses of both planets are randomized between 0.25 Mj 
and 4 Mj. Similarly, for the 5:3 resonance, the masses 
are randomized between 0.25 Mj and 5 Mj, and the 
2:1 between 1 Mj and 12 Mj. The initial phase of the 
inner planet is set to zero, while the outer planet's is 
randomized between and 27r. 

To determine the initial position of the outer planet, 
we set the planet's semi-major axis so that it is just out- 
side the particular resonance by a small fixed amount 
(« 0.1 AU). We then apply a frictional force of the 
form — av, where a is a constant and v is the veloc- 
ity of the outer planet. It is applied for a short amount 
of time to allow the outer planet to migrate into the 
resonance location with the inner planet. The time 
and strength have been adjusted so no significant cou- 
pled migration occurs. Once the planets begin to mi- 
grate together in resonance, the change in their semi- 
major axis beco mes correlated with the cha nge in their 
eccentricity (see iMurrav fc Dermott I I2OOO0 . Since no 
significant coupled migration occurs, the eccentricity of 
each planet never grows above values of 0.2 from migra- 
tion/eccentricity excitation alone. 

To monitor the resonance in question, we use the res- 
onance variable defined by 

V = jlK + 32K + js^o + ji^i (3) 

(jMurrav fc Dermott |[2000l) . Here A and uj are the mean 
longitude and longitude of pericenter, respectively. The 
subscripts denote either the inner or outer planet. Since 
our integrations are two-dimensional, the terms includ- 
ing the longitude of the ascending node (f2) are zero. For 
a p : q resonance, the ordered set {ji, J2, J3, J4} is equal 
to {p, -q,-{p-q), 0} or {p, -q, 0,~{p~q)}. When two 
planets are in a given resonance with pericenter align- 
ment, the resonance variable is bounded with a libration 
amplitude < 2tt. Additionally, planets can still be locked 
in resonance although ip continues to oscillate between 
and 27r. Such configurations require oscillations of both 
the semi-major axis and eccentricit y to account for the 
circu lating longitude of pericenters ( Murrav fc Dermott I 
[2OOOI) . 

A run is defined as stable if the system evolves for 
10^ yr and maintains the original configuration of semi- 
major axes without appreciable changes in any orbital 
elements (eccentricity, etc.). In these runs, the planets 
must lock into a mean motion resonance so that cj) be- 
comes bounded with a finite libration amplitude or circu- 
lates in a periodic manner. A run is defined as unstable 



if the system goes unstable before or after the planets 
lock into mean motion resonance. 

3.2. Results 

We present the results for the resonance stability tests, 
analyzing both the region of stability for these particu- 
lar resonances as well as the final configuration of the 
systems that went unstable. We first focus on each reso- 
nance separately, and will consider broader results in the 
final section. 

3.2.1. The 2:1 Resonance 

We determined analytically from Equation [2] that the 
2:1 resonance is stable up to Mi + Mj « 10 Mj. Figure 
[T] presents a stability map, where the filled circles rep- 
resent a stable run, and the (red) open circles represent 
runs where instability occurred before the planets could 
lock into a 2:1 MMR. The dashed line in the figure is the 
analytical boundary between stable and unstable regions 
determined by Equation[2l It includes all three terms and 
hence is not symmetric in the masses. It's slope therefore 
changes depending on whether the inner or outer planet 
is the more massive one (recall the definition of in 
Equation [2|) . The upper left region of the plot shows 
an obvious region of stability not predicted by the sim- 
ple analytic criterion. That is, the additional stability 
typically involves a larger outer planet, and includes up 
to two 11 AIj planets remaining stable. However, there 
is another obvious region of instability not predicted by 
Equation [2] for an inner planet roughly twice the mass 
of the outer planet. In these configurations, the planets 
went unstable due to the inner planet exciting the ec- 
centricity of the outer planet. This resulted in the outer 
planet undergoing a close encounter with the inner planet 
and typically was ejected from the system. 

Over 90% of the unstable runs resulted in the less mas- 
sive planet being ejected from the system, leaving the 
remaining planet with eccentricities that typically range 
from 0.3 all the way to 0.97. The latter of these would 
be expected to undergo strong tidal forces at periastron 
and result in a less eccentric orbit with a smaller semi- 
major axis (see Section HJ. The time scales to undergo 
instability are typically very short. Most of the runs that 
go unstable before entering resonance do so within 10** 
yr. Once the planets entered resonance to a point where 
(j) become bound and pericenter alignment occurred, the 
planets enjoyed a phase protection that prevented any 
further close encounters. 

3.2.2. The 3:2 Resonance 

Equation[2]suggests that the 3:2 resonance is stable up 
to Ali + Mj w 2 Mj. This region of stability is plotted in 
Figure [2l which uses the same conventions as Figure [TJ 
Additionally, the (red) crosses represent runs where the 
planets showed signs of being in a 3:2 MMR but even- 
tually went unstable. The region of stability for the 3:2 
resonance is more symmetric compared to the 2:1 reso- 
nance, with the region of stability extending to over twice 
the analytically determined values. The 3:2 resonance 
typically does not include any configurations where one 
of the planets is more than 3 Mj. 

The majority of the instabilities result in a collision, 
leaving the remaining planet in a position between the 
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TABLE 1 

Mass Limits for particular resonances as determined by Equation \2\ Here p denotes 

THE ratio of the OUTER AND INNER PLANET'S PERIOD, To/Ti. LIMITS ARE GIVEN DEPENDING ON 
WHETHER ONLY ONE (fIRST-ORDER) , TWO, OR ALL THREE TERMS ARE USED IN EQUATION [2] 



Resonance 


P 


Ml + M2 


Ml + M2 


Mx + Ah 


I'irst-Order / 






{Mj, first-order) 


{Mj, second-order) 


(Mj, third order) 


Third-Order 


3:2 


1.50 


2.26 


1.54 


1.77 


1.272 


5:3 


1.66 


5.05 


3.13 


3.72 


1.356 


2:1 


2.00 


15.33 


8.09 


10.10 


1.517 


5:2 


2.50 


45.15 


19.66 


25.86 


1.745 


3:1 


3.00 


95.29 


35.57 


48.58 


1.961 







initial semimajor axes of the two original planets with a 
very low (< 0.1) eccentricity. In the several cases where 
the instability did not immediately result in a collision 
and the two planets strongly interacted for an extended 
period of time, the simulation resulted in either a colli- 
sion or ejection, both leaving the remaining planet with 
a higher final eccentricity. A low number of runs (< 5%) 
resulted in the planets being thrown into a configuration 
that did not go unstable for the remainder of the of the 
integration. In these cases, the ratio of periods, T2/T1 
remained in the interval [5.0,9.0]. 

There exists a few cases where planets begin to enter a 
3:2 MMR but eventually go unstable. These instabilities 
occur in less than 10^ yr. An example is shown in Figure 
[21 The top panel plots one of the resonance variables: 



ip = pXo - qXt. - (p - q)uJo, 



(4) 



where the middle panel plots difference in pericenters of 
the two planets, and the bottom panel plots the eccen- 
tricity of the inner and outer planet (black and red lines, 
respectively). We see that the planets show sign of en- 
tering resonance around 5000 yr, where both ip and Alu 
show periodic variation from that point onward. How- 
ever, shortly after 1.1 • 10"* yr, the two planets leave the 
resonance and go unstable due to strong gravitational 
perturbations. This particular run resulted in a collision. 

3.2.3. The 5:3 Resonance 

Finally, Figure [4] shows the stability map for the 
second-order resonance^, the 5:3. Compared to the 2:1 
and 3:2 MMRs, the region of stability is less than that 
predicted by analytic criteria with only three stable runs 
existing beyond the region of analytically determined sta- 
bility. These three runs, however, are quite close to the 
border between the two analytically determined regions. 
Additionally, a number of unstable runs exist deep within 
the region predicted to be stable, showing that configura- 
tions in a 5:3 resonance are very sensitive to their initial 
conditions. The majority of the simulations result in an 
ejection, leaving the final planet with an eccentricity in 
the interval [0.3,0.8]. Collisions result in outcomes sim- 
ilar to that of the 3:2 resonance: a merger product with 
low eccentricity somewhere between the semimajor axes 
of the two original planets. 

The above three resonances are potentially very im- 
portant for any pair of planets that form in reasonably 
close proximity to each other. Planets can be captured 
into these resonances without having appreciable initial 
eccentricity, making them highly likely candidates for 

All resonances can be written in the form (a + b)/a, where b 
denotes the order of the resonance. 
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Fig. 1. — StabiUty map for planets initially on circular orbits 
in the 2:1 resonance. A filed circle shows a run where the plan- 
ets remained in the resonance for lO'^ yr, and an open (red) circle 
shows runs that went unstable before entering a 2:1 MMR. Plan- 
ets were determined whether to be in resonance by monitoring the 
resonance variable defined by Equation \3\ The analytical bound- 
ary between stability and instability determined by Equation [2] is 
shown by the dashed line, with the analytical region of stable lying 
below the line, and unstable otherwise. 

planets within the protoplanetary disk, where disk ef- 
fects are believed to initially keep planets on circular or- 
bits. These resonance impose a strict geometrical con- 
figuration on the orbits of the planets in and near these 
resonances. These periodic interactions between the two 
planets can lead to instabilities, resulting in ejections or 
collisions. Figure [5] plots the cumulative distribution of 
the remaining planet's eccentricity after an ejection or 
collision occurred due to an instability that formed dur- 
ing the integration. 

Figure O compares the eccentricity distribution with 
that of observed extrasolar planets (solid red line). Ejec- 
tions (dotted black line) produce a similar shape dis- 
tribution, but underproduce lower eccentricity planets. 
Collisions (dashed black line) produce the majority of 
planets with eccentricities below 0.1. In order to com- 
pare these various processes with observed eccentricities, 
we randomly select 50 outcomes (which can be both colli- 
sions or ejections) from each of the three resonances. We 
then plot that distribution (solid black line) in Figure [SI 
and find it matches observed eccentricities better than 
collisions or ejections alone. For further comparison, we 
include the results of non- resonant planet-planet sc atter- 
ing (red dotted hue) from lChatteriee et al. I ()2007l ). 
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Fig. 2. — Stability map for planets initially on circular orbits in 
the 3:2 resonance. A filed circle shows a run where the planets re- 
mained in the resonance for 10^ yr, and an open (red) circle or cross 
shows runs that went unstable before or after entering a 3:2 MMR, 
respectively. Planets were determined whether to be in resonance 
by monitoring the resonance variable defined by Equation \3\ The 
analytical boundary between stability and instability determined 
by Equation [2] is shown by the dashed line, with the analytical 
region of stable lying below the line, and unstable otherwise. 
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Fig. 3. — Evolution of two planets that enter a 3:2 MMR but 
eventually become unstable and result in a collision. The top panel 
plots one of the resonance variables, defined by Equation |3] while 
the middle panel plots the pericenter difference and the bottom 
panel plots the eccentricity of the inner and outer planet (black and 
red lines, respectively). The two planets enter resonance around 
5000 yr, but strong gravitational perturbations break the resonance 
around 1.1 ■ 10* yr. 
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Fig. 4. — Stability map for planets initially on circular orbits in 
the 5:3 resonance. A filed circle shows a run where the planets re- 
mained in the resonance for lO'^ yr, and an open (red) circle or cross 
shows runs that went unstable before or after entering a 5:3 MMR, 
respectively. Planets were determined whether to be in resonance 
by monitoring the resonance variable defined by Equation |3] The 
analytical boundary between stability and instability determined 
by Equation [2] is shown by the dashed line, with the analytical 
region of stable lying below the line, and unstable otherwise. 
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Fig. 5. — Cumulative distribution showing the final eccentricity 
of the remaining planet after instability. It includes the data from 
all three resonances (the 3:2, 5:3, and 2:1) without the presence of 
disk forces. Those due to collisions and ejections are shown by the 
dotted black and dashed black lines, respectively. The solid red 
line is the eccentricity distribution of observed extrasolar planets. 
It omits systems which are believed to be in resonance or are close 
enough to their parent star to undergo tidal interactions. The solid 
black line is a randomly selected combination of 50 outcomes from 
each of the three resonances. This is used to graph the overall dis- 
tribution from all possible outcomes. It is compared with the final 
eccentric ity distribution du e to no n-resonant planet-planet scatter- 
ing from lChatteriee et al. I 1120071) . shown by the dotted red line. 
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4. PLANET-PLANET-DISK INTERACTIONS 

In this section we model planets interacting with a 
protoplanetary disk. The previous section considered 
the stability of three resonances without the effects of 
planet-disk interactions. We now study the stability of 
these resonances while the planets undergo migration due 
to angular momentum exchange between the planet and 
the surrounding disk material. Doing so will incite fur- 
ther eccentricity growth for the two planets, which can 
greatly affect the stability of the system. 

We model the protoplanetary disk as a pure gas disk 
with an inner hole that extends all the way to the central 
star. A self-consistent treatment of planet-disk interac- 
tion would require us to also model the feedback of the 
planet on the disk, which in turn would necessitate incor- 
porating hydrodynamics into our simulation. To simplify 
matters, we impose the "closing in" of the inner hole in- 
dependently of the two planets' posit ions, but at a rate 
consistent with results from t heory ([Ward I Il997f) an d 
full hydrodynamic simulations (jOgilvie fc Lubowl 2003f ) . 
Therefore, the disc extends only to the location of the 
outer planet, and the inner planet essentially sits in a 
"gap," and does not undergo migration unless it is due 
to resonance interactio ns. Hydrodynamic simulations 
(|Thommes et al. 1 l2008t ) have shown that systems with 
multiple gap-opening planets tend to quickly remove all 
of the material inside the outermost planet, which puts 
the configuration in exactly the state we have assumed. 

The "radius" of the disk's inner edge is set to migrate 
linearly inward at a rate of 1 AU every 10^ yr. A Fermi 
function is then used to model the edge of the disk, al- 
lowing the friction effects to smoothly fade away. That 
is, given the "radius" of the disk edge, r^, the accelera- 
tion due to friction applied to each planet is modeled by 
the following equation: 



af 



-F{a,rd)-[ — 



(vr)r 



where F(a,rd) is the Fermi function 



F{a,rd) = 1 



1 



l + exp[(a-rd)/(0.05rd)] 



(5) 



(6) 



for a given semi-major axis a. The "width" of the Fermi 
function (O.OSr^) is selected so only the outer planet feels 
non-negligible force from the protoplanetary disk unless 
the two planets are going to unavoidably collide.^ As the 
inner disk edge contracts, the planet adjusts its position 
relative to the disk edge until it reaches an equilibrium 
location. At this location, the torque it experiences is 
equal to that needed to make it migrate inward at the 
same rate as the disk edge. 

The two terms in Equation [5] model the migration and 
eccentricity damping, respectively. The migration time 
scale, tm, is set to 6000 yr, consistent w ith the type I rate 
for a body with mass of order 10^ Mq ()Ward Ill997f) : the 
exact value is unimportant as long as it is short com- 
pared to the time scale of the disk edge contraction. The 
eccentricity damping time scale, te, is set equal to either 

More specifically, it is chosen so the inner planet is always four 
"widths" away from the radius of the disk. It is easy to show that, 
at this distance, the inner planet will always feel less than 2% of 
the full magnitude of the frictional effects. 



a tenth the value of t™, as suggested by previous ana- 
lytical aiid_jiumericaJ_^^ 

disks (iGoldreich fc Tremainc "19 8GHTrimng et aU llQQa 
Nelso n et all HoOO: Papaloizou c t al. Il2001h or 10^° yr, 
so to model the disk without eccentricity damping. 

As planets locked in resonance migrate together, their 
change in eccentricity is related to the mass ratio and 
rate of change in the sem i-major axis of the planets 
(jMurrav fc Dermott Il2000l ). For planets that migrate to 
less than a tenth of their original semi-major axis, appre- 
ciable eccentricity growth occurs. Therefore, our simu- 
lations take into account energy lost due to the close-in 
tidal interactions between eccentric planets and the cen- 
tral star. As a planet approaches from large distances 
with large eccentricities, the relative distance is short at 
periastron compared to the rest of the orbit. During this 
brief period, angular momentum and energy is exchanged 
between the star and planet. We adopt the app roxima- 
tion developed bv lPapaloizou fc Terquem] ()2001[ ). where 
the tidal interaction acceleration is written as 



iGrUpRl 



0.6m 



where i?, is the stellar radius and Rp {— j^/2GM^,) is 
the distance of closest approach for a parabolic orbit with 
angular momentum j. This model is an approximation of 
a small perturbation limit for a non-rotating star and the 
planet being on a parabolic (or highly eccentric) orbit. 
It neglects the effects of tides acting on the planet it- 
self. Since these interactions include a term proportional 
to the velocity of the planet, changes in the semi-major 
axis occurs rapidly and the planet "falls" into the central 
star on short time scales. Therefore, we do not include 
relativistic effects in our simulations since the planets 
do not spend considerable time at distances were these 
effects would be appreciable. 

4.1. Initial Conditions 

The initial semi-major axis of the inner planet is ran- 
domized between 5 and 6 AU, and the mass of both plan- 
ets is fixed at 1 Mj. The outer planet's semi-major axis is 
then randomized between the minimal distance at which 
the planets are stable according to Equation [5] (which is 
very near at the 3:2 resonance for two Mj planets) and 
slightly outside the 2:1 resonance, that is. 



Oi + 2.4 (/Zi -I- /Lto) 



1/3, 



< Oo < 22/3a, 



(8) 



where e is roughly 0.1 Oi and /j, is defined as in Equation 
[3 The initial location of the "radius" of the disk edge is 
chosen so that the outer planet lies approximately four 
"widths" away from the radius, which allows the outer 
planet to travel along the outer edge of the disk as it 
migrates inward without any sudden changes in the ex- 
ternal disk forces. 

Two separate sets {N « 500) of integrations are done, 
one with te = O.ltm and the other with t^ — 10^° yr 
(hereafter referred to as t^ = oo). 

4.2. Results 

As is to be expected given the initial condi- 
tions described above — particularly the low initial 
eccentricities — all planets are captured into one of three 
resonances: 2:1, 5:3, and 3:2. An important parame- 
ter is the amount of eccentricity damping due to the 
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disk, which affects the overaU stabihty of the planets 
in these resonances. We find that the 2:1 resonance is 
the most stable; all planets captured into it remain in 
the 2:1 until they have migrated close enough to the star 
for tidal interactions to pull the planets apart (dragging 
the inner planet into the star). The 3:2 is similar for 
the case where damping is present and eccentricities are 
kept small. When eccentricity damping is not present, 
the 3:2 resonance goes unstable at eccentricities larger 
than 0.3. The 5:3 resonance goes unstable for very low 
eccentricities and produces the largest variety of possible 
outcomes with and without eccentricity damping. 

Each run has four possible results. The planets can 
either (1) capture into a particular resonance and mi- 
grate together until the inner planet is close enough to 
tidally interact with the star and the resonance is bro- 
ken, resulting in the inner planet colliding with the star. 
In this case, we classify the run as "Stable;" (2) cap- 
ture into a particular resonance and undergo instability, 
so that the two planets collide with each other. The 
few cases (« 1%) where instability results in one of the 
planets colliding with the star are also included here; 
(3) undergo instability once in a resonance and result in 
one of the planets being ejected; (4) after migrating in a 
particular resonance, undergo instability but ultimately 
end up in a different resonance. This occurs because one 
planet is thrown into the disc and is driven back inward 
towards the disc edge, recapturing the inner planet into 
a new resonance. Once in this resonance, the planets ei- 
ther quickly become unstable or remain in the resonance 
for remainder of the integration. The cases that remain 
in this new resonance are labeled "Scattered and Recap- 
tured." If they quickly become unstable and result in a 
collision or ejection, they are classified in the appropriate 
category (2) or (3) above. 

We examine each case separately, with and without 
damping. The results are given in Table [2l 

4.2.1. Overall Results with te/tm = 0.1 

We first present the case where eccentricity damping 
is present, which is set to be reasonably representative of 
planet-disk interactions (see Section 2] above) . The 2:1 
and 3:2 resonances are always stable. However, the 5:3 
resonance goes unstable at very low eccentricities. Over 
half of the 5:3 resonances end up with a collision, re- 
sulting in a single planet with low eccentricity around 
6 AU (recall that we stop integrations as soon as there 
is a collision or ejection). Since the planets have equal 
mass, about half of these cases eject the outer planet, 
and half the inner. The 5:3 resonance, however, was able 
to produce planets in a wide variety of resonances that 
included the 2:1, 3:2, 3:1, 5:3, 4:1, 6:1, and even a few 
cases of planets in a 10:1 and a possible 17:1 (where both 
planets held eccentricities above 0.95). These resonances 
were monitored by calculating the appropriate resonance 
variable defined by Equation [31 These latter resonances 
were produced by one of the planets being scattered out- 
ward due to instability. With this planet now possessing 
a significant eccentricity, it was able to be captured into a 
higher-order resonance as it resumed inward migration. 
Both planets typically reached very high eccentricities 
during their coupled migration in these new resonances. 
Examples of a stable system and a scattered and recap- 
tured system are given in Figures [6] and [T] 
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Fig. 6. — The dynamical evolution of two planets in a 2:1 res- 
onance with the presence of disk and tidal forces, letting te/tm = 
0.1. Top Panel: Plot of the semi-major axes of the outer and inner 
planet, shown by the red and black lines, respectively. The blue 
lines show the location of the disk's inner edge, with the outer blue 
linos being the width of the edge, and the middle line being the 
disk's inner radius. Middle Panel: A plot of the absolute value of 
the difference in longitude of pericenters (black), and the two 2:1 
resonance variables defined by Equation |3] (red and magenta lines). 
The resonance breaks apart at the very end due to the presence of 
tidal forces acting on the inner planet, which pull the planets apart. 
Bottom Panel: The eccentricity of the inner and outer planet, us- 
ing the same color convention as the top panel. The decrease in 
eccentricities near the end of the evolution is due to the presence 
of tidal forces from the star, modeled by Equation [7] 




2x10= 4x10= 6x10= 8x10= 10= 

t (yr) 

Fig. 7. — The dynamical evolution of two planets that go unstable 
in a 5:3 resonance and get trapped into a 6:1 resonance. This occurs 
with the presence of disk and tidal forces, assuming te/tm = 0.1. 
Top Panel: A plot of the semi-major axes of the outer and inner 
planet, shown by the red and black lines, respectively. The blue 
lines plot the location of the disk's inner edge, with the outer blue 
lines being the width of the edge, and the middle line being the 
disk's inner radius. Middle Panel: A plot of the absolute value of 
the difference in longitude of pericenters (black), and the two 6:1 
resonance variables defined by Equation |3] (red and magenta lines). 
The resonance breaks apart at the very end due to the presence of 
tidal forces acting on the inner planet, which pull the planets out 
of the resonance. Bottom Panel: The eccentricity of the inner and 
outer planet, using the same color convention as the top panel. 
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4.2.2. Overall Results Without Eccentricity Damping 

We now repeat the simulations with eccentricity damp- 
ing completely removed, allowing us to see how the out- 
comes look in the limit of maximal eccentricity excita- 
tion. They are also presented in Table [2] for easy com- 
parison. 

There is no significant change in the 2:1 and 5:3 res- 
onances. The biggest change comes with the 3:2 reso- 
nance, which now rarely survives the migration towards 
the star. The majority of the 3:2 cases result in a colli- 
sion. Also, ejections are more common for the 3:2 reso- 
nance than for the 5:3. The 3:2 also has a few cases where 
the planets are scattered and recaptured into another 
resonance, the possible resonances comprising a smaller 
subset of those mentioned above. 

4.2.3. The 5:3 Resonance 

The 5:3 resonance is worth discussing in more detail. 
For the case where eccentricity damping is present, we 
tally the distribution of eccentricities of the remaining 
planet after a collision or ejection has occurred. This 
are shown in Figure [H Additionally, Table [3] gives the 
various resonances (with relative capture probabilities) 
produced by instabilities with the 5:3 resonance. 

Collisions leave the remaining merged planet with very 
low eccentricities, although there is a larger range in 
their semi-major axes compared to our earlier integra- 
tions (See Figure[5]). Roughly 30% of the planets are left 
closer than 5 AU, the minimum initial semi-major axis 
of the inner planet. The rest are found between 5 and 7 
AU. 

Ejections tend to leave the remaining planet at less 
than 3 AU with eccentricities above 0.6, and a large per- 
centage of these planets having eccentricities between 0.6 
and 0.8. This is similar to the results from where the 
eccentricity distribution approximately spanned the in- 
terval [0.5^0.8]. 

Similar to the procedure used in 331 we also compare 
an overall distribution, but this time focusing entirely on 
the 5:3 resonance. Table [3] shows that in the case with 
eccentricity damping, the ratio of collisions to ejections 
is 6:1. Therefore, we randomly select 25 integrations re- 
sulting in an ejection, and 150 integrations resulting in a 
collision and plot the distribution in Figure O It is clear 
that the 5:3 resonance alone cannot accurately reproduce 
the observed eccentricity distribution in extrasolar plan- 
etary systems below eccentricities of 0.5. 

5. SUMMARY AND DISCUSSION 

In this paper, we have examined the stability of partic- 
ular first and second-order resonances with and without 
disk interactions. We have shown that the two-planet 
stability criterion (computed to third-order) is invalid for 
these resonances (see Figures [H [21 and|4|). The 2:1 res- 
onance has regions of stability that extend up to mass 
ranges comparable to brown dwarfs, while the 3:2 res- 
onance has a region of stability with combined planet 
masses up to twice the above criterion. The 5:3 reso- 
nance has a reduced region of stability compared to an- 
alytic criteria. 

Next, we have examined the results of instability. As 
shown in Figure [51 combining an equal number of un- 
stable outcomes originating in each of the above three 




e 



Fig. 8. — Cumulative distribution of the remaining planet's ec- 
centricity after a collision or ejection occurred due to an instability 
in the 5:3 resonance (black dashed and dotted-lines, respectively). 
The solid red line shows the distribution of observed extrasolar 
planetary systems. The selection method is the same as in Figure 
[5] The black solid line is a randomly selected distribution includ- 
ing both collisions and ejections, which included the results of 25 
integrations that resulted in an ejection and 150 integrations that 
resulted in a collision. The fraction taken from each was deter- 
mined by their relative percentages of occurrence, which are given 
in Table [2] 



TABLE 2 

Relative percentages for the final possible outcomes 
given that two 1 mj planets are caught into a particular 
resonance and undergo disk-planet interactions. cases 
with and without constant eccentricity damping are 
presented. a run is "stable" if the planets remain in that 
resonance throughout the entire integration, resulting in 
one of the planets migrating into the star. cases where 
AN "Ejection" or "Collision" occurred are noted 

APPROPRIATELY. If THE PLANETS UNDERGO INSTABILITY AND END 

up locking in another resonance for the majority of the 
integration, it is labeled as a "scattered and 
Recaptured." 



Resonance 


Final Outcome 


te/tm = 0.1 


= OO 


2:1 


Stable 


100% 


100% 




Collision 


0% 


0% 




Ejection 


0% 


0% 




Scattered and Recaptured 


0% 


0% 


3:2 


Stable 


100% 


3% 




Collision 


0% 


68% 




Ejection 


0% 


27% 




Scattered and Recaptured 


0% 


2% 


5:3 


Stable 


< 1% 


0% 




Collision 


54% 


54% 




Ejection 


9% 


12% 




Scattered and Recaptured 


36% 


34% 







resonances (with appropriate relative fractions of ejec- 
tions and collisions) yields an eccentricity distribution 
which fits, at least, as well as previous studies of dy- 
na mical instability in mul ti-planet systems (in particu- 
larlChatteriee et al. II2007I for the three-planet case, and 
I Juric fc Tremaine I l2007l for larger ensembles). This is 
an intriguing result which raises the possibility of a link 
between these MMRs and exoplanet eccentricities. How- 
ever, significant caution is warranted in interpreting it. 
Our "outcomes" represent snapshots at the time of in- 
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TABLE 3 

Possible resonances that two 1 Mj planets were thrown 

into after they locked into a 5:3 resonance and 
ultimately went unstable due to eccentricity excitation 
and disk-planet interactions. the results are presented 
for the case where damping was present. resonances are 

included if coupled migration results after the two 
planets fall into these particular resonances such that 
the ratio of periods remained fixed for the migration. 
Additionally, we monitor that the corresponding 
resonance variables, defined by equation [3] librate with 

amplitude less THAN 2lT OR UNDERGO SYSTEMATIC PERIODIC 

LiBRATiON. See Figure [7] for an example. The overall 

RELATIVE PERCENTAGES OF CAPTURE FOR EACH RESONANCE IS 
GIVEN. 



Order 


Resonance 


Relative % 


First-Order 


2:1 


21.3% 




3:2 


1.7% 


Second-Order 


3:1 


37.1% 




5:3 




Third-Order 


4:1 


4.9% 




5:2 


13.5% 


Other 


5:1 


1.1% 




6:1 


3.1% 




7:2 


< 1% 




7:1 


5.1% 




8:1 


2.3% 




9:1 


2.8% 




10:1 


2.8% 




11:1 


1.7% 




12:1 


1.4% 




14:1 


< 1% 




17:1 


< 1% 







stability, with no attempt made to model the subsequent 
effect of the disk on the remaining single planet. Also, 
the way in which we assemble resonances in i 33.1l is rather 
artificial, at least in comparison to convergent migration 
of SI 

With the addition of migration in 21 plus the eccen- 
tricity damping expected from disk interactions, we find 
the 2:1 resonance to be the most stable, with 100% of the 
runs remaining stable even without eccentricity damping 
(see Table [2]). Additionally, the 3:2 resonance remains 
completely stable with our adopted damping prescrip- 
tion. A second-order resonance, the 5:3, has a high cap- 
ture efficiency but goes unstable at very low eccentrici- 
ties due to close encounters. This means that in order for 
planets to become locked into the 3:2 resonance, their pri- 
mordial period separation must generally lie between 5:3 
and 3:2. In this study, we selected the initial semi-major 
axes of the planets to coincide wit h the location where it 
is believed planets form (see, e.g.. lKokubo fc Ida1l2002l : 
iThommes et al. ]|2003[ ). rather than allowing the plan- 
ets to migrate from further out in the disk. Under this 
assumption that planets form between 5 and 10 AU, if 
planets were able to form at any location with equal prob- 
ability, the probability that a planet forms between the 
3:2 and 5:3 resonances of a planet located at 5 AU is 



7r( [(5/3)^3,(3/2)^/^] -5) 
7r-52 



0.009. 



(9) 



That is, the probability of forming such a configuration 
is about 1%. 

To date, we have not discovered any planets in a closer 
than 2:1 resonance. The lack of 5:3 planets is readily 
accounted for by this study: Since this resonance is un- 
stable at low eccentricity, it is unstable to survive all 



the way to a mature planetary system. The fact that 
3:2 planets are also not seen, notwithstanding the res- 
onance's robustness, may well be telling us something 
fundamental about how planets form: It may simply be 
that neighboring giant planets are seldom born with pri- 
mordial period separations of less than 5:3. 

The planets in HD 12661 are believed to be in a 6:1 
resonance, and there is much debate as to how they could 
have become locked in such a high-order resonance. This 
study suggests a mechanism. We have shown that an 
instability in one of the lower-order resonances, which 
needs little to no initial eccentricity for planets to be- 
come trapped in it, can result in the planets scattering 
apart, to then be brought together again and re-locked 
in a more distant, higher-order resonance (See Table [3]). 
Further work is necessary to include the adjustment of 
the disk edge after the planets have been locked into a 
more distant MMR. However, migration of the outward- 
scattered planet would likely be somewhat more gentle in 
a self-consistent disk — generally we would expect an an- 
nulus of gas to initially exist between the two planets — in 
which case resonant re-capture would actually tend to be 
more likely than in our current simple implementation. 

Aside from the physical findings, this study highlights 
an important point on orders of accuracy. Since the res- 
onance stability tests relied on two planets of unequal 
mass, we chose to use up to the third order term in Equa- 
tion [21 because this was the first term that was not sym- 
metric in the masses (that is, a mapping Mi Mo would 
not be an identity mapping). Furthermore, as shown by 
Table 1, additional terms can change the result up to 
twice the first-order value. If we were to truncate Equa- 
tion [21 to just the first term, we would have concluded 
that the equation gives a better representation of the re- 
gion of stability for the resonant case than it actually 
does. 

This paper represents a first step in our study of reso- 
nances in protoplanetary discs. Here we have constructed 
systems where resonance encounters were impossible to 
avoid: requiring that we assumed the planets were copla- 
nar and existed near these resonances without going un- 
stable. A natural question this paper does not address is 
whether the dynamics that might occur due to relative 
inclinations affects the ove rall stability of the resonant 
systems. It has been shown (jThommes 'fc Lissauer II2003D 
that it is possible for planets starting in a 2:1 resonance 
to later enter an inclination resonance, which gener- 
ates large mutual inclinations even in an initially almost 
coplanar system. It is not clear a priori whether this 
would promote dynamical instability or protect against 
it. Additionally, we have modeled the eccentricity damp- 
ing time scale as a constant pr oportional to the migra- 
tion time scale. Recent stud ies (jOgilvie fc Lubow"1l2003t 
iMoorhead fc Adams I l2007t l have shown that damping 
may not be uniform and depends on the local properties 
of the disk as well as the overall geometry of the planets. 
In particular, de/dt may even change sign depending on 
the value of the eccentricity. Therefore, it is important to 
better understand the nature of eccentricity growth and 
damping in order to understand the long-term stability 
of these resonances in the protoplanetary disk. 

Our improved distribution shown in Figure [5] was de- 
rived by using an equal number of cases from each of 
the three resonances studied. However, Table [2] suggests 
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that we should find more planets in the 2:1 resonance 
compared to the 3:2 resonance during the disk's lifetime 
when we assume that all the planets migrated from a re- 
gion farther out than the 2:1 resonance. In the case of 
gravitational instability, it may be assumed that plan- 
ets could form in a more equally distributed manner be- 
fore undergoing migration and locking into these reso- 
nances. Complete three dimensional integrations would 
allow us to gain an accurate representation of the dis- 
tribution of resonance captures due to migration alone. 
Further stability tests will involve allowing one planets 
to migrate from greater distances, where near-resonance 
interactions cannot affect the overall stability of the sys- 
tem until the planets (given they remain stable) can mi- 
grate close enough together to encounter a particular res- 
onance. 

Additionally, gap-forming planets on the edge of a disk 
in general ought to still be accreting some mass across 
their gap, and therefore mass growth will be taken into 
consideration. In this case, planets could lock into reso- 
nance under conditions we have shown to be stable, and 
through mass accretion enter a regime that analytic and 
numerical criteria determine to be unstable. A compari- 
son of whether this region of stability is similar to what 
we have found here for constant-mass planets would give 
us further important insights into the possible range of 
planetary system configurations. 

We have gained a clearer picture of the range of possi- 
ble dynamics within a protoplanetary disk. In doing so, 
we have seen how resonances can strongly affect plane- 
tary systems, long before the surrounding disk material 
has dissipated. This suggests that what happens before 
the disk disappears is key in determining a planetary 



system's ultimate dynamical architecture. One possible 
outcome is a system which, after the disk is gone, is left 
with planets in resonance. In order for this to happen, 
the planets must enter these resonances relatively late in 
the disk's lifetime, so that they ultimately survive the 
coupled migration. After the disk has dissipated, some 
resonances may eventually become unstable. Another 
possibility we have demonstrated is that resonances can 
be broken by dynamical instability that occurs while the 
disk is still present. Our results allow us to make a con- 
nection between the planet formation process and obser- 
vations of mature, resonant planetary systems, none of 
which have yet been observed in a closer than 2:1 reso- 
nance: For planets that form further apart than a period 
ratio of 2:1 and migrate convergently, the very stable 2:1 
constitutes a formidable barrier. Between 2:1 and 5:3, 
the initial outcome is likely to be capture into 5:3, fol- 
lowed by a one-way trip to instability. The 3:2 resonance 
is quite stable (unless eccentricity damping by the disk is 
absent), so provided that an appreciable number of plan- 
ets form with a primordial period separation of less than 
the 5:3, we would expect to find some planets in a 3:2 res- 
onance. Thus, if future observations continue to reveal 
no such systems, then this may suggest a lower limit to 
orbital separations with which neighboring planets first 
form. 

This work was supported by NSF Grant AST-0507727 
at Northwestern University. ET would like to acknowl- 
edge the support of the Spitzer Theoretical Research Pro- 
gram, and NSERC of Canada. 



REFERENCES 



Barnes, R. & Greenberg, R. 2007, ApJ, 665, L67 
Bodenheimer, P., Hubickyj, O., & Lissauer, J. J. 2000, Icarus, 143, 
2 

Boss, A. P. 2000, ApJ, 536, 101 

Cameron, A. G. W. 1978, Moon Planets, 18, 5 

Chatter jee, S., Ford, E. B., & Rasio, F. A. 2007, Accepted to ApJ 

(arXiv:astro-ph/0703166v2) 
Fischer, D. Marcy, G. W., Butler, R. P., and Vogt, S. S., Henry, 

G. W., Pourbaix, D., Walp. B., Misch, A. A., & Wright, J. T. 

2003, ApJ, 586, 1394 
Ford, E. B., Havlickova, M., & Rasio, F. A. 2001, Icarus, 150, 303 
Ford, E. B., & Rasio, F. A. 2007, Accepted to ApJ 

(arXiv:astro-ph/0703163v2) 
Gladman, B. 1993, Icarus, 106, 247 
Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 
Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 
Goldreich, P., & Sari, R. 2003, ApJ, 585, 1024 

Juric, M., & Tremaine, S. 2007, Accepted to ApJ 

(arXiv:astro-ph/0703160v2) 
Kley, W., Peitz, J., & Bryden, G. 2004, A&A, 414, 735 
Kokubo, E., & Ida, S. 2002, ApJ, 581, 666 
Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 59 

Lee, M. H., Butler, R. P., Fischer, D. A., Marcy, G. W., & Vogt, 

S. S. 2006, ApJ, 641, 1178 
Lin, D. N. C, & Papaloizou, J. C. B. 1979, MNRAS, 186, 799 
Lin, D. N. C, & Papaloizou, J. C. B. 1986, ApJ, 309, 846 
Lin, D. N. C, &; Papaloizou, J. C. B. 1993, in Protostars and 

Planets III ed. E. H. Levy &; J. I. Lunine (Tucson: University of 

Arizona Press) p. 749 
Lissauer, J. J. 1993, ARA&A, 31, 129 

Marcy, G. W., & Butler, R. P. 1995, 187th AAS Meeting BAAS, 
27, 1379 

Marcy, G. W., & Butler, R. P. 1998, ARA&A, 36, 57 

Marcy, G. W., Butler, R. P., Fischer, D., Vogt, S. S., Lissauer, J. J., 

& Rivera, E. J. 2001, ApJ, 556, 296 
Mayor, M., & Queloz, D. 1995, Nature, 378, 355 



Moorhead, A. V., & Adams, F. C. 2007, Icarus, 708 

Murray, C. D., & Dermott, S. F. 2000, Solar System Dynamics 

(Cambridge: Cambridge University Press) 
Nelson, R. P., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, 

MNRAS, 318, 18 
Ogilvie, G. I., & Lubow, S. H. 2003, ApJ, 587, 398 
Papaloizou, J. C. B., & Lin, D. N. C. 1984, ApJ, 285, 818 
Papaloizou, J C. B., & Terquem, C. 2001, MNRAS, 325, 221 
Papaloizou, J. C. B., Nelson, R. P., & Masset, F. 2001, A&A, 366, 

263 

Papaloizou, J. C. B. 2003, Celest. Mech. Dyn. Astro., 87, 53 
Papaloizou, J. C. B., Nelson, R. P., Kley, W., Masset, F. S., 

Artymowicz, P. 2007, in Protostars and Planets V ed. B. 

Reipurth, D. Jewitt, & K. Keil (Tucson, Univ. of Arizona Press) 
Press, W. H., Teukolsky, S. A., VetterUng, W. T., & Flannery, B. P. 

1992, Numerical Recipes in C. The art of scientific computing 

(2nd ed.; Cambridge University Press) 
Rasio, F. A., & Ford, E. B. 1996, Science, 274, 54 
Snellgrove, M., Papaloizou, J. C. B., & Nelson, R. P. 2001, A&A, 

374, 1092 

Thommes, E. W., & Lissauer, J. J. 2003, ApJ, 597, 566 
Thommes, E. W., Duncan, M. J., & Levison, H. F. 2003, Icarus, 
161, 431 

Thommes, E. W., Matsumura, S., & Rasio, F. A. 2008, Science, 
808 

TriUing, D. E., Bnez, W., Guillot, T., Lunine, J. I., Hubbard, W. 

B., & Burrows, A. 1998, ApJ, 500, 902 
Udry, S., Fischer, D., & Queloz, D. 2007, in Protostar and Planets 

V ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson, Univ. of 

Arizona Press) 
Ward, W. R. 1986, Icarus, 67, 164 
Ward, W. R. 1997, Icarus, 126, 261 

Weidenschilling, S. J., & Marzari, F. 1996, MNRAS, 180, 57 



